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PREFACE 


Titanium  matrix  composites  reinforced  with  SIC  fibers  have  been  the  subject  of  extensive 
research  over  the  past  decade  because  of  interest  in  advanced  aerospace  applications  that 
require  light  weight,  high  stillhess,  and  high  temperature  capabilities.  Programs  such  as  the 
National  Aerospace  Plane  (NASP)  and  the  Air  Force  Integrated  High  Performance  Turbine 
Engine  Technology  (IHPTET),  in  particular,  have  been  evaluating  the  properties,  mechanical 
behavior,  and  failure  mechanisms  of  this  class  of  materials  under  typical  anticipated  use 
conditions.  While  failure  modes  depend  on  the  specific  mechanical  and  thermal  loading 
conditions,  fiber  failure  has  been  identified  as  either  a  direct  mode  of  failure  or  as  a  dominant 
mode  of  failure  after  some  other  mode  such  as  matrix  cracking  has  produced  an  increase  in  fiber 
stresses. 

Direct  fmlure  of  fibers  which  act  as  a  bundle  embedded  in  a  ductile  matrix  material  is 
considered  to  be  the  primary  cause  of  composite  material  failure  under  creep,  low  cycle  fatigue, 
and  in-phase  thermo-mechanical  fatigue  (IP  TMF).  The  process  of  breakage  of  the  fiber  bundle, 
however,  is  not  very  well  understood.  Results  of  IP  TMF  tests  on  composites  with  SiC  fibers 
having  different  matrix  materials,  and  tested  at  different  maximum  temperatures,  show  apparently 
different  modes  of  progressive  failure  after  a  fiber  has  broken.  Timetal  matrix  composites  at  650 
C  maximum  temperature  show  debonding  of  the  fiber/matrix  interface  in  the  region  adjacent  to 
the  fiber  crack.  Ti-6-l-4V  composites  at  427  C,  on  the  other  hand,  show  an  apparent  propensity 
for  matrix  cracking  adjacent  to  the  fiber  break.  To  evaluate  the  potential  for  fi^rming  an  interface 
crack,  and  to  shed  light  on  the  mechanics  of  crack  formation  and  interfacial  toughness,  the  present 
study  was  undertaken.  The  analysis  refrains  from  making  many  of  the  simplifying  assumptions 
about  stress  distributions  which  are  commonly  utilized  in  analyses  of  this  type.  Instead,  a  rigorous 
and  systematic  elastic  analysis  is  conducted  using  only  the  boundary  conditions  as  an"assumption" 
in  addressing  the  problem.  Moreover,  the  general  problem  is  broken  down  into  key  problems 
which  are  examined  individually  and  will  finally  be  combined  to  provide  answers  to  the  general 
problem. 

It  was  found  convenient  to  divide  the  work  into  two  different  parts. 

In  parti,  we  examine  the  load  transfer  characteristics  of  a  material  system  after  a  fiber 

breaks. 


In  part  II,  we  examine  the  residual  stresses  due  to  curing  and  thermal  stresses  due  to 
differences  between  the  thermal  expansion  coefficients  of  the  matrix  and  fiber. 

In  part  III,  we  are  developing  a  sinc-function  based  numerical  method  that  is  geared 
specifically  for  solving  3D  problems  that  posses  stress  Angularities,  in  certain  regions  of  the  body. 
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by 

£.  S.  Folias 
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ABSTRACT 


In  this  study,  the  load  transfer  characteristics  of  a  broken  fiber  are  investigated. 

The  problem  consists  of  a  cylindrical  fiber  that  is  embedded  into  a  matrix  material.  The  fiber  axis 
is  assumed  to  coincide  with  the  z-axis  and  a  crack  is  assumed  to  be  present  on  the  plane  z  =  0  and 
for  r  <  a.  Far  away  fi-om  the  crack,  the  fiber  is  subjected  to  uniform  external  load  of  o„. 

Moreover,  adjacent  to  the  crack  and  along  the  interface  ( see  fig.  3  ),  the  matrix  and  fiber 
surfaces  are  assumed  to  slide  along  the  interface  length  -c  <  z  <  c,  in  the  presence  of  a  variable 

fiiction  T/z  =  .  On  the  other  hand,  perfect  bonding  is  assumed  to  prevail  all  along  the 

remaining  interface,  i.e.  for  |  z  |  >  c. 

For  the  solution  of  the  problem,  we  utilize  a  Fourier  Integral  Transform  whereby  we 
reduce  the  problem  to  that  of  the  solution  of  a  singular  integral  equation  along  the  interface  path 
<  <  z  <  c.  The  solution  of  this  integral  equation  then  allows  the  determination  of  the  displacement 
and  stress  fields.  Two  areas  of  special  interest  immediately  come  to  mind,  ( i )  the  neighborhood 
adjacent  to  the  point  |  z  |  =  c  and  r  =  a,  and  ( ii )  the  neighborhood  adjacent  to  the  point  z  =  0  and 
r  =  a. 


The  analysis  reveals  that  the  load  transfer  characteristics  depend  heavily  on  the  material 
properties  of  the  system. 


1.  Introduction 


It  is  well  recognized  that  realistic  modeling  of  deformation  and  damage  accumulation  in 
the  fatigue  of  metal-matrix  composites  depends  heavily  on  the  mechanisms  which  govern  the 
damage.  Experimental  evidence  of  the  accumulation  of  strain  in  a  SiC-fiber/titanium-matrix 
composite,  in  conjunction  with  fractographic  examination  of  the  samples,  point  to  fiber  breakage 
as  the  dominant  mode  of  failure  during  thermomechanical  fatigue  tests.  During  this  physical 
process  individual  fibers  break  periodically,  at  random  locations,  and  throughout  the  composite. 
Considerable  work  on  the  modeling  of  broken  fibers  in  composites  has  already  been  conducted 
over  a  number  of  years.  In  the  literature  one  may  find  the  work  of  Rosen  (1964)  using  a  shear-  lag 
model.  The  stress  concentration  present  in  the  immediate  vicinity  of  a  broken  fiber  was  first 
analyzed  by  Hedgepeth  (1961)  and  has  subsequently  been  addressed  by  many  others  who  have 
utilized  other  simpler  models.  Recently,  Penado  and  Folias  (1989)  investigated  the  stress 
concentrations  in  a  fiber  based  on  three  dimensional  considerations.  Moreover,  the  3D  edge 
effects,  when  a  fiber  meets  a  free  surface,  have  been  investigated  by  Folias  (1989).  A  closer 
inspection,  however,  of  the  local  neighborhood  of  a  fractured  fiber  (  see  Fig.  1)  reveals  that  the 
problem  is  much  more  difficult  that  it  was  originally  thought  of  This  is  because  one  must  also 
account  for  the  vertical  interface  cracks  that  develop  along  the  sides  of  the  fiber.  One  may  also 
recall  the  fact  that  composites,  in  general,  are  attractive  in  using  them  for  practical  applications 
because  the  fibers  are  designed  to  carry  most  of  the  applied  load.  Consequently,  when  a  fiber 
breaks,  at  least  locally,  the  load  must  then  be  transferred  somehow  to  the  other  portion  of  the 
broken  fiber  via  the  matrix.  This  suggests,  therefore,  that  the  adjacent  matrix  must  now  carry  a 
much  higher  stress  load.  In  reality,  however,  one  may  conjecture  that  ( i)  there  is  some 
redistribution  of  the  load  to  the  adjacent  unbroken  fibers,  and  that  (ii)  there  is  perhaps  an  increase 
in  the  load  which  the  adjacent  matrix  must  now  carry.  Recently,  Nicholas  and  Ahmad  (1994), 
through  the  analysis  of  a  relatively  simple  model,  were  able  to  provide  some  qualitative  answers 
regarding  this  complex  phenomenon.  On  the  other  hand,  a  more  sophisticated  stress  analysis  of 
the  local  problem,  coupled  with  the  results  of  a  sporadic  sequence  of  fiber  breaks  along  the 
composite  material  ^stem,  can  provide  important  information  and  guidance  to  material  designers 
for  the  development  of  future  material  systems. 

2.  Formulation  of  the  Problem 

Consider  the  equilibrium  of  a  material  system  that  occupies  the  space  |x|<oo,|y|<oo  , 
I  z  I  <  00  and  contains  a  cylindrical  fiber  that  is  embedded  into  a  matrix  space  (  see  Fig.  1  ).  The 
fiber  axis  coincides  with  the  z-axes  and  the  fiber  radius  is  r  =  a .  Both  fiber  and  matrix  materials 
are  assumed  to  be  homogeneous,  isotropic  and  linearly  elastic.  At  the  interface  ( r  =  a  )  perfect 
bonding  is  assumed  to  prevail. 

In  this  model,  we  envision  the  failure  to  take  place  in  three  different  stages.  During  the 
first  stage,  a  plane  crack  is  formed  in  the  middle  of  the  fiber  region,  i.e.  at  r  <  a  and  on  the  plane 
z  =  0,  and  advances  until  it  reaches  the  fiber/matrix  interface  boundary.  This  problem  has  been 
studied  by  others  and  recently  more  rigorously  by  Pagano  et  al  (1995 ).  Subsequently,  during  the 
second  stage,  a  certain  length  along  the  fiber/matrix  interface,  above  and  bellow  the  plane  z  =  0, 
has  debonded  and  the  matrix  is  now  allowed  to  slip  in  the  presence  of  a  non-constant  fiictioiuil 
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force.  Finally,  during  the  third  stage,  the  fiber  crack  pups  open  and  a  finite  jump  in  the  fiber 
displacement  is  formed.  The  present  analysis  deals  with  the  events  of  the  second  stage.  We  believe 
that  the  analysis  of  this  stage  will  provide  important  information  related  to  the  failure  process  and 
may  reveal  the  sufficient  conditions'  required  for  the  subsequent  advancement  of  the  failure  to 
that  of  stage  three.  Putting  it  another  way,  it  is  hoped  that  it  will  reveal  the '  sufficient  conditions ' 
required  to  suppress  any  further  development  of  the  failure. 

At  a  certain  instant  in  time,  the  fiber  breaks  in  two  peace^s  shown  in  Figs2  and  3.  The 
separation  distance  A©  ,  denotes  the  maximum  displacement  between  the  fiber  faces  after  the 
system  has  been  loaded.  It  is  also  assumed  that,  subsequent  to  the  fiber  breakage,  splitting 
between  fiber  and  matrix  along  the  interface  and  up  to  a  distance  -c  <  z  <  c  is  simultaneously 
being  formed  in  the  presence  of  a  non-uniform  frictional  force.  As  to  loading,  far  away  fi-om  the 
location  of  the  break,  the  fiber  carries  a  uniform  stress  a©  parallel  to  the  direction  of  the  fiber 
axis.  Similarly,  the  matrix  carries  a  uniform  stress  ai  also  in  the  direction  of  the  fiber  axis. 

In  the  absence  of  body  forces,  the  differential  equation  governing  the  stress  potential 
functions  and  is  given  by 

V^O=0  (1) 

where  is  the  biharmonic  operator.  The  displacement  and  stress  fields  are  then  given  in  terms 
of  the  potential  functions  as  ; 

(i)  displacement  field : 


(2) 

2Gw=2(1-v)V2<D-0 

(3) 

(  a  )  stress  field : 

(4) 

<T^.  =  {2-v)V=f-0 

(5) 

(6) 

where  V  is  Poisson's  ratio  and  G  the  shear  modulus. 

As  to  boundary  conditions,  we  require  that : 
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(i)  matrix  region : 

at  z  =  0  :  =  0  ( 7 ) 

atz  =  0:  tS’^  =  0  (8) 

(ii)  fiber  region: 

atz  =  0:  tS  =  0  (9) 

at  z  =  0 :  Jo  2Ttrc^^dr=  0  ( lo ) 

atz  =  0:  vi^^(r,0)  =  Aoi/a^-r^  (ii) 


It  may  be  noted  here  that  the  latter  represents  a  crack  that  has  the  shape  of  a  very  sharp 
ellipse.  In  this  study,  the  primary  emphasis  was  to  examine  the  effect  that  a  non-uniform  frictional 
force  has  on  the  mechanism  of  failure  and  to  gain  some  further  insight  on  the  evolution  of  the 
events  of  this  complex  crack  formation.  In  a  follow  up  paper,  which  the  authors  are  presently 
working  on,  the  boundary  conditions  have  been  relaxed  by  allowing  now  the  fiber  displacement  w 
to  possess  a  finite  jump. 

(Hi)  interface  region : 


at  r  =  a : 

u?  = 

u?”> 

;0<z<oo 

(12) 

at  r  =  a : 

Cn  = 

eg” 

;0<z<oo 

(13) 

at  r  =  a : 

Vy(0  = 

;0<z<oo 

(14) 

at  r  =  a : 

Trz  = 

Trz 

;  c<z<oo 

(15) 

at  r  =  a : 

TfZ  = 

\lGrr 

;  0  <  z  <  c 

(16) 

Perhaps  it  may  be  appropriate  here  to  note  that  one  of  the  difficulties  presented  by 
boundary  condition  (  16  )  is  that  the  shear  stress  Xrz  is  an  odd  fiinction  of  z  while  the  normal 
stress  is  an  even  function  of  z .  At  first  glance  this  appears  to  be  O.K.,  however  the 
mathematical  satisfaction  of  such  type  of  a  boundary  condition  requires  some  very  special 
attention. 


(iv)  far  away  from  the  plane  z  =  0 : 

Czz  =  CTo 


(17) 
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(18) 


where  Oi  is  related  to  Go  via  the  relation  : 

lilTl|zl-w(Ezz  ~  =  0  (  19  ) 

Moreover,  in  order  to  complete  the  formulation  of  the  problem  we  must  require  fh^t : 
as  r  &  z-^oo  :  all  complementary  displacement  and  stresses  must  vanish 
and  the  continuity  condition  : 

at  r  =  0  ;  all  the  complementary  displacements  and  stresses  for  the  fiber  must  be 
bounded 


3.  Method  of  Solution 

\^^thout  going  into  the  mathematical  details,  we  have  constructed  the  following  stress 
potential  functions  that  exhibit  the  proper  behavior  at  infinity. 

(si)}s\n(sz)ds 

±  Z  {  C?  +  D?an  Izl  }exp(-an \z\)JoianO  ( 20 ) 

n=0 

+aV>z^  +  b^^zf^ 

where  the  upper  sign  refers  to  z  >  0  and  the  lower  sign  refers  to  z  <  0 , 

den  =  2GfVm  -  GmVf-  GfVf(1  -  2Vm)  -Gf+Gm-  2VfVmGm  (  21  ) 


~  ~  ~  ^Gfiym  +  Vf)  +  4V/nV/Gf 

+Gf-Gni}  (22) 

“  2(den)  ~  2Vn7)  +  V/nG/n(1  —  2vr)}  (  23  ) 

0^'”>  =  Jo  {/\('”>Ko(s/)  +  B^"’\sl)K^{s^)}sm(sz)ds 

+a('”)z^  +  b^'^^zr^  ( 24 ) 
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with 


=  “eSo  HVinVf(Gf  -  Gm)  -  4VfGf  +  2Gm(Vf  -  Vm) 

+G,-6™}^  (25) 

-  V/Gf + 2v„Vf(Gf-  G„)}  ^  ( 26 ) 

The  constants  are  to  be  determined  from  the  boundary 

conditions  of  the  problem. 

By  construction,  boundary  conditions  (  7  )-(  8  )  are  automatically  satisfied.  Next,  to  satisfy 
boundary  condition  (  9  ),  we  let 

C?  =  2vrD?  (  27 ) 

Similarly,  to  satisfy  boundary  condition  (  10  )  we  choose  ,  n  =  0,1,2,... ,  to  be  the  roots  of 
the  equation 


*^l(otn3)  — 0  (28) 

in  view  of  which  eq.  (  10  )  now  reduces  to 

-v,)S®+A«]/i(sa)+saB»/o(sa)}(/s=-fa.  (29) 

We  will  suppress,  at  this  time,  the  satisfaction  of  the  above  condition  but  will  return  to  it  at  a  later 
time. 


Next  from  eq  (  1 1  ),  we  have 

V9<')(r.  0)  =  -2^  £  alDfjo(a„n  =  Ao  Va^ 

n=n 


in  view  of  which  the  coefficient  D^n  may  now  be  determined  as 

=  n  =  1.2.3,.,. 

(ana)  2 


Gf  jg(a„a)  3  - 


dnd 


■^iD2>  =  |A„a 


(30) 


(31a) 


(31b) 
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Thus,  the  crack  opening  in  the  fiber  region  becomes  a  veiy  sharp  ellipse.  Perhaps  it  may  be 
appropriate  here  to  note  that  the  constant  A©  has  now  been  related,  through  eqs  (  29  )  and  (31), 
to  the  applied  load  Go  .  Moreover,  the  matrix  load  stress  CTi  is  related  to  the  fiber  load  stress 
CTo  through  the  eq.  (  19  ).  More  specifically, 

_  f  Gm  h^^myrGf+Gm+2vmVfGr2GmV/Vm+G/Vm+GmVm] 

^  '  Gf  l2Gt\'m-GmVrG/v/+2vmy/Gf-Gf+Gm-2v/VmGm]  ‘  (32) 

Next,  the  boundary  conditions  at  the  interface  r  =  a 

l/P-ur  =  0  (33) 

aff  -  =  0  ( 34 ) 

become  respectively, 

Jo  s2  {(sa)/Ce(sa)e('”)  +  K^  (sa)A('”)  +  G(sa)/o(sa)B(') 

+G/i  (sa)A^^  }  cos(sz)c/s  =  0  (35) 

and 

s2{[/,(sa)  -  (sa)/„(sa)]>4<'>  +  [-(1  -  2v,)(sa)/„(sa) 

-(sa)^l^(sa)]B^^  +  [Ki(sa)  +  (sa}Ko(sa)]A^^ 

+[-(1  -  2vm)(sa)Ko(sa) + {sa)^K^  (sa)]S^'”)  }cos(sz)c/s  = 

=  -  f  aa^ { 1  - ttn Izl Jo(ana)exp(-an Izl)  ( 36 ) 

/M) 

Similarly,  the  boundary  condition' 

Vl^O  _  _  Q  ^  27  ) 

becomes 

Jo  ( lo{sa)A^^  +  [4(1  -  Vf)lo(sa)  +  (sa)/i  (sa)]B^^ 

-GKo(sa)/\<'”>  +  [4(1  -  v„)GKo(sa)  -  G(sa)/<1  (sa)]S<'”) }  ( 38 ) 

xsin(sz)cfs=  S  (2(1  -  Vf)  +  anlzl}a5D?Jo(ana)exp(-anlzI) 

/}=0 

'  See  remark  on  top  of  page  3. 
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where  we  have  adapted  the  definition 


G  = 


Gm 


In  order  to  facilitate  our  subsequent  discussion,  we  define 


e  =  (sa) 


lo(sa). 

/i(sa)> 


=  (sa) 


Ko(sa) 

Ki(sa) 


and 


^(0  = 


A(f)  . 
h(say 


B(0  = 


B(0 

/i(sa) 


Ko(sa)  ’ 


D(m)  _  Jai. 

Ki(sa) 


(39) 


(  40  a,b  ) 


(41a,b) 

(42a,b) 


In  view  of  eqs  (  40  a,b  ),  (  41  a,b  ),  (  42  a,b ),  equations  (  35  )  -  (  36  )  and  (  38  ),  by  Fourier 
inversion,  become  respectively : 


A(f)  +  =  0  ( 43 ) 

(1  -  e)A(^  +  [-(1  -  2vf)e  -  (sa)2]B(0  +  (i  + 

+[-(1  -  2v;;,)e,  +  (sa)2]B('”)  =  |f  J*  sin(sOcf^  ( 44 ) 

and 

e-A®  +  [4(1  -  v,)e + (sa)2]B®  -  GeiAf"* 

+[4G(1  -  v,)e,  -  G(sa)2]S("’>  =  If  |^{2(1  -  vf)U 

-|f}Sin(S4)C(^  (45) 

where  we  have  made  use  of  the  definition 

U(4)=  Sa|D?exp(-a„l5l)Jo(a),a)  .  (46) 

n=0 

A 

The  system  of  eqs.  ( 43  )-  (  46  )  may  now  be  solved  in  terms  of  one  of  the  unknowns,  say  B^'"^  . 
More  specifically. 
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(47) 


Am  =  +  Ifi  ||sin(s|)c/5 + |f^  Usin(s|)d| 

D(D  —  [(1-G)(sa)^-(1-6)(1-2vm)®H-2G(1-Vni)®f]  DCm') 

°  -  2(1-vr)®  ^  2(W  ° 


+ll?  Jo  ^sin(s^)cf^ 


(48) 


=  -ea(^  - 

where  for  simplicity  we  have  adopted  the  following  definitions  : 

P,  =  1  -  G + e,  4  Ge + -  G^i  ^  $) 


P,  =  -(sa)^  4  (1  -  2v.)e,  -  G(1  -  e)e,  - 

^J(1-G)(sa)2-(1-G)(1-2vn,)®f+-2G(1-Vm)®f] 

X  ^ 

o  [(1-2vr)®+«(1-®)+{sa)2] 

^'3  = - 5 - 

Finally,  the  last  boundary  condition  becomes 


Ilmr-^axS 


limr-»a  Xrz  ^  I 

limr-j-aCucrSf^) 


C<Z<oo 

0<z<c 


(49) 

(50) 

(51) 

(52) 

(53a) 

(53b) 


where 


tS  =  Jo  { [2(1  -Vf)+e]B^^  +  A^^} exp(-s(a -  0)sin(sz)ofs  ( 54 ) 

=  -J^ +  [-2(1  - vn,) + ef]B('”>}exp(-s(r- a))sjn(sz)c/s  ( 55 ) 
cf?’'  =  -l/^s2{[1 4e,]A<"’)4Kl  -2v„)e,4(sa)2]e<'”)}  (56) 

xexp(-s(r-  a))cos(sz)cfs 

It  remains,  therefore,  for  us  to  satisfy  the  last  boundary  condition  (  53  ),  as  well  as  eq. 

(  29  ).  This  will  ultimately  tie  the  remaining  constant  in  terms  of  the  applied  load  of  the  material 
system. 
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4.  The  Integral  Equation 


In  order  to  solve  the  dual  integral  equation  (  53  ),  we  first  cast  it  into  the  form  of  a  Cauchy 
singular  integral  equation  of  the  first  kind.  For  this  reason,  we  let  at  r  =  a 


JJ*  s^(sa)B('”>sin(sz)cfs  =  \\f(z)  ;  0  <  z  <  c 

(57) 

and 

s^(sa)B('”>cos(sz)c/s  =  (j)(z)  ;  0  <  z  <  c 

(58) 

where  Vj/  represents  an  odd  function  of  z  and  (j)  represents  an  even  function  of  z .  Furthermore, 
the  reader  may  be  reminded  that,  by  construction,  equation  (  S3a  )  is  automatically  satisfied  along 
the  interface  segment  C  <  Z  <  oo  .  Thus,  it  remains  for  us  to  satisfy  equation  (  53b )  along  the 
interface  segment  0  <  Z  <  C.  Without  going  into  the  mathematical  details,  equation  ( 53  )  may 
now  be  written  as,  after  some  straight  forward  manipulations. 

aC"^f«’'d5  =  /i:2)  ;  0<z<co 

(59) 

where  we  have  adopted  the  following  definitions  : 

(60) 

with 

f  f  4b3(6+n+2(ivr)  246® 

(z2+62)3  ^ 

(61) 

,  -4+2vrG  62(-16+8vr4G-12n)  24m6^  ^ 

(z2+6*2)3  (z2+62)4  /  • 

(62) 

In  writing  the  function  f  (z),  we  have  made  use  of  the  following  approximation 

(63) 

with 


^=^(3^;  Ao  = -0.2089^^1^';:  b  =  0.0700a  (64abcd) 

4G(1-vO 
“  ”(3-4vffG)  • 

The  general  solution  of  the  above  integral  equation  is  given  by 
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K 


(65) 


where  k  now  represents  an  arbitrary  constant  and  where  the  integrals  are  to  be  evaluated  in  the 
Cauchy  Principal  Value  sense.  Without  going  into  the  tedious  mathematical  details,  the  integrals 
in  equation  (  65  )  may  now  be  evaluated  and  the  following  result  is  obtained  after  simplification. 


\l/(Z)-M(t)(Z) 


2  e  Aq  a  f  cfeP  c^o^z  -v 

I(zz0r 


(66) 


where  the  coefficients  Cle\do^  ,  k  —  0,  1, 2, 3, 4  are  long  expressions  involving  the  material 
constants  of  the  system.  Thus,  in  view  of  equations  (  65  ),  (  57  )  and  (  58  ),  one  may  now 

determine  explicitly  the  coeflRcient  .  In  particular. 


s\sa)Bf’"^  =  -^A,U 

/c—0 


4^  C  rJliSC) 
"1  1 


kc^  J2(sc)  - 
(c2+62)  (SC)  +---J 


I  1  r«/o(SC)  kc^  4i($C) 

^  (c2+t»2)»fL  1  +  (SC) 


(67) 


Moreover,  in  view  of  equations  (  47  )  -  (  49 )  and  (  67  ),  the  remaining  coefficients  may  now  be 
determined  and  the  displacement  and  stress  fields  can  be  recovered  explicitly. 

5.  The  Stress  Fields 


As  it  was  prewously  noted,  it  is  now  relatively  easy  to  recover  the  displacement  and  the 
stress  fields  everywhere  within  the  fiber,  and  within  the  matrix,  region  of  the  material  system.  Two 
areas  of  special  interest  are  worthy  of  examination,  ( i )  the  neighboiiiood  of  the  point  |  z  |  =  c  and 
r  =  a  (  see  Fig.  4  ),  and  ( ii )  the  neighborhood  of  the  point  |  z  |  =  0  and  r  =  a.  In  this  report,  we 
examine  the  stress  field  of  the  former.  Suppressing  the  long  and  tedious  algebraic  manipulations, 
one  can  show  the  stresses  to  be: 


stress  field  at  I  z  I  =  c  and  r  =  a : 

tS  =  coBo^l^  {cos(f)  +  MSin(|)}  ( 68 ) 

+croBi  yi”  sin(<|)){cos(f )  +  nsin(f )}  +  0(8°) 

=  croBeyj{cos(|)  -  IlSin(|)}  (  69  ) 

+<ioBi  yf  sin(<j)){cos(|)  -  Iisin(f )}  +  0(8°) 
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(70) 


=  -o<,6o^{ncos(j)+ sin(|)} 

-ooBi  yj  sin((|>){^cos(f ) + sin(f )}  +  0(8°) 
eg'  =  -eg’  +  0(8'’)  ( 71 ) 


where 


p  3.032f)  A  1 

Do-  c  ^Oool  (3-4V/+6)  ^ 


Bi  = 


(1-G)  p 
4(1-vr)°o 


(72b) 


Finally,  returning  to  eq  (  29  ),  we  see  that  the  quantity  Ao  may  now  be  related  to  the 
applied  load  Go  ■ 

6.  Conclusions 

Without  going  into  the  numerical  detdls,  it  would  be  of  practical  interest  to  examine  how 

well  the  actual  physical  boundary  condition  CTk  =  0  is  satisfied  on  the  plane  z  =  0.  For  this 
reason,  we  assume 

^=10,n  =  0.5,f=5,v,  =  0.25  (73) 

and  upon  defining 

one  finds  that 

GfSAo  +  csor  =  0  ,  (  75  ) 


or  upon  solving  for  Ao 

=  =  0.629^  .  (76) 

In  view  of  the  present  analysis,  the  numerical  result  for  the  function  S  is  given  in  Fig.  5. 
The  reader  wUl  notice  that  its  profile  is  a  straight  line  up  to  approximately  r/a  =0.65  where  by  it 
begins  to  deviate.  This,  however,  was  to  be  expected  for  it  is  the  result  of  our  approximation  of 
the  function  U  (see  equation  (63))  with  a  single  term.  Moreover,  it  may  be  noted  that  the  function 
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S  does  not  exhibit  any  oscillations.  While  this  approximation  does  eflFect  the  boundary  condition 

aS  =  0  in  the  neighborhood  of  the  point  r  =  a  and  z  =  0,  it  can  be  shown  that  it  has  no  effect 
the  singular  term  of  the  stress  fields  in  the  vicinity  of  the  crack  tip  |  z  |  =  c.  In  conclusion,  the 
actual  boundary  condition  is  seen  to  be  satisfied  a  lot  better  than  expected. 

Returning  next  to  the  stress  fields,  we  note  that  the  usual  singular  stress  behavior, 

which  is  characteristic  to  crack  problems,  still  prevails  in  the  vicinity  of  the  crack  tip,  i.e.  at  the 
interface  point  |  z  |  =  c.  Furthermore,  the  unexpected  condition 

aT  +  cP  =  0{z°)  ill) 

also  holds,  suggesting,  therefor,  that  an  increase  in  the  stress  a  a  is  followed  by  a  proportional 
decrease  in  the  stress  . 

As  a  practical  matter,  the  analysis  will  now  be  used  to  derive  a  criterion  with  which  to 
predict  the  debonded  interface  length  2c,  as  well  as  the  sufficient  conditions  required  to  suppress 
any  further  development  of  the  interface  crack  to  that  of  stage  three  (  please  see  discussion  on  top 
of  page  2  ).  Moreover,  as  it  was  previously  noted,  we  are  presently  examining  the  third  stage 
where  the  fiber  crack  has  now  pupped  and  where  the  fiber  displacement  w  has  been  allowed  to 
posses  a  finite  jump.  The  results  will  be  reported  in  a  follow-up  paper  in  which  the  effects  of  a 
thermal  loading  have  also  been  included . 
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Fig.  1  Geometrical  configuration  for  perfect  bonding  and  no  friction  at  the  interface. 


matrix 


Fig.  2  Geometrical  configuration  for  a  cracked  fiber  with  no  friction  at  the  interface. 


fiber 


U2 


1/ 


Fig.  3  Geometrical  configuration  of  broken  fiber  with  friction  at  the  interface. 
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fiber  /  matrix  interface 


region  of  fiiction 


I 
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Fig.  4  Local  coordinate  system  at  the  point  z  -  c. 
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j{azz  -  Go}  z=odr 


±=10,h  =  0.5,|  =  5,v,=  0.25 


r/a 


Fig.  5  Examination  of  the  actual  boundary  condition  Gzz  on  the  plane  z  -  0. 
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PREDICTING  CRACK  INITIATION  IN  COMPOSITE  MATERIAL 
SYSTEMS  DUE  TO  A  THERMAL  EXPANSION  MISMATCH 


by 

E.  S.  Folias 
Michael  Hohn 


ABSTRACT 


Residual  stresses  due  to  curing  and  thermal  stresses  due  to  differences  between  the 
thermal  expansion  coefficients  of  the  matrix  and  fiber  may  have  a  major  effect  on  the 
microstresses  within  a  composite  material  system  and  must  be  added  to  the  stresses  induced  by 
the  external  mechanical  loads.  Such  microstresses  are  often  sufficient  to  produce  microcracking 
even  in  the  absence  of  external  mechanical  loads,  example  during  the  cooling  process. 

In  this  report  a  few  selected  results  are  presented  for  a  material  system  consisting  of  SIC-6  - 
cylindrical  fibers  which  are  periodically  embedded  into  a  plate  matrix  consisting  of  beta21 
material.  The  results  are  based  on  a  linear  elastic  micromechanics  model  which  provides  the  stress 
profiles  due  to  (i)  a  uniform  load  perpendicular  to  the  direction  of  the  fibers  and  (ii)  due  to  a 
thermal  expansion  mismatch.  In  this  analysis,  perfect  bonding  between  the  fiber  and  the  matrix  is 
assumed  to  prevail.  For  this  case  the  analysis  shows  that  at  a  free  edge  and  for  a  lateral  load,  there 
exists  a  weak  stress  singularity  which  increases  as  the  temperature  increases.  Selected  stress 
profiles  are  given  for  the  above  two  loads.  Moreover,  the  application  of  a  fracture  criterion  shows 
that  no  failure  is  likely  to  take  place  for  a  cooling  temperature  AT  of  900°  C.  Thus,  the  growth  of 
any  pre-existing  microcracks  will  be  suppressed. 


1.  INTRODUCTION. 


Residual  stresses  due  to  curing  and  thermal  stresses  due  to  differences  between  the 
thermal  expansion  coefficients  of  the  matrix  and  fiber  may  have  a  major  effect  on  the 
micro-stresses  within  a  composite  material  system  and  must  be  added  to  the  stresses  induced  by 
the  external  mechanical  loads.  Such  micro-stresses  are  often  sufficient  to  produce  micro-cracking 
even  in  the  absence  of  external  loads,  example  during  the  cooling  process.  Furthermore,  if  the 
material  system  is  thermally  fatigued,  these  residual  stresses  may  cause  some  of  the  existing 
micro-cracks  to  grow  and  coalesce  and  thus  form  the  presence  of  larger  cracks. 

Thus,  if  rational  designs  in  the  use  of  fiber-reinforced  metal  matrix  composites  are  to  be 
made,  their  performance  under  static,  dynamic,  and  thermally  fatigued  loads  need  to  be 
predictable.  The  first  step  towards  this  goal  is  the  realization  that  the  ultimate  failure,  as  well  as 
many  other  aspects  of  the  composite  behavior,  are  the  result  of  growth  and  accumulation  of 
microdamage  to  the  fibers,  matrix  and  their  interfaces.  Thus,  it  appears  that  any  generally 
successful  model  of  performance  and  failure  must  incorporate  the  effects  of  this  damage  in  some 
way.  This  certainly  represents  a  challenge.  In  this  paper,  we  address  the  form  of  such  damage  due 
to  the  residual  stresses  developed  as  a  result  of  the  thermal  expansion  mismatch  between  the 
fibers  and  the  matrix. 

In  this  work,  a  systematic,  3D,  micromechanics  approach  is  used  in  which  the  fibers  of  a 
composite  material  system  are  modeled  as  cylindrical  inclusions  that  are  embedded  into  a  matrix 
plate.  The  analytical  model  is  then  used  to  predict,  the  residual  stresses  due  to  a  thermal  expansion 
mismatch,  e  g.  during  a  cooling  process.  Moreover,  the  model  provides  a  better  understanding  of 
how  the  residual  stresses  are  developed  and  how  they  can  be  controlled  particularly  in  relation  to 
ceramics  where  there  is  no  ductility  to  accommodate  any  plastic  deformation. 

The  analysis  reveals  the  dependence  of  the  residual  stress  field  on  the  fiber  volume 
fraction  ratio,  identifies  the  critical  locations  where  a  crack  is  most  likely  to  initiate  and 
subsequently  propagate,  recovers  the  interface  shear  stress  profile  and  provides  important 
information  and  guidance  to  material  designers  for  the  pre-selection  of  fiber  and  matrix  materials 
in  order  to  alleviate  some  of  the  residual  stresses.  It  may  be  noted  that  the  theoretical  model  is 
applicable  to  ceramic  and  metal/matrix  composite  systems. 

2.  MATHEMATICAL  MODEL. 

Consider  in  infinite  plate  matrix  which  consists  of  material  Beta21,  see  fig.  1 .  The  matrix 
plate  is  assumed  to  extend  to  infinity  both  in  the  x-  and  y-  directions.  In  the  z-direction,  the  matrix 
plate  is  assumed  to  have  a  finite  dimension,  2h,  in  order  to  capture  any  possible  3D  effects  that 
may  be  present.  A  uniform  and  square  periodic  set  of  cylindrical,  SIC-6,  fibers  are  einbedded  into 
the  matrix  plate  in  the  directions  of  both  x  and  y.  Two  types  of  loads  are  being  considered  :  (i)  a 
uniform  transverse  load  perpendicular  to  the  direction  of  the  fibers  and  along  the  y-direction 
and  (ii)  a  uniform  temperature  load  AT  (cooling )  that  is  applied  throu^out  the  material  system. 
Both  fibers  and  matrix  are  assumed  to  be  homogeneous  and  linear  elastic. 
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The  governing  equations,  are  the  well  known  Navier's  equations  coupled  with  the  Energy 


Balance  equation.  More  specifically. 

i-2vat  1  V  «  a,. 

(1) 

1  ^  +  V2v  -0 

l-2vdy  ^  l-2v  ^ 

(2) 

1  y2^  2{1+v)^_q 

l-2v  az  ^  l-2v  dz 

(3) 

S7^T=0 

(4) 

As  to  boundary  conditions,  (i)  the  appropriate  stresses  are  required  to  vanish  at  the  free 
edge,  (ii)  perfect  bonding  is  assumed  to  prevail  at  the  fiber  /  matrix  interface,  (iii)  finally  the 
boundary  conditions  within  the  cell  configuration  are  required  to  be  satisfied.  Once  the 
displacement  field  has  been  completely  determined,  the  stresses  can  be  obtained  by  using  the 
stress-strain  relations: 

Oy  =  X,5yEy  +  2Gzij  -  a(3X  +  2G)(T-  T a)5ij. 

(5) 

3.  THE  3D  DTSPLACEMENT  FIELD, 


Without  going  into  the  mathematical  details,  the  3D  displacement  field  has  been  derived  by 
Folias  (  1976)  and  can  be  expressed  as  : 


^  ^  -  ll/'iCPvr)  +m/2(pvz)} 


mj-2^^  dx 


(6) 


+ 1  ^  cos(an/j)cos(a„r)  +lf 


d^/{> 


rt=Q 


^  S  ^{2(W;  -  lyi(Pvz)  +  m/2(Pvz)} 


(7) 


V=1 


00  arrC^)  ,  .  S/Wy-l  Jj)  I  2 

-I  ^cos(a„/j)cos(a„z)  +  -;;;^/3  +^2  a,2 

n=l  ^ 


dif 


=  -^  E  ^{-2(Wy  -  ll/'lCPvZ)  +  w/2(3vZ)}  -  (  ») 

V=:l 


where. 


/i(Pvz)  =  cos(Pv^)cos(Pvz) 

/2(Pvz)  =  (pv/j)sin(Pv/j)cos(PvZ)  -(pvz)cos(pv/j)sin(Pvz) 


(9) 

(10) 


3 


(11) 


+ 

(^+^-p;)«S>=o.  (12) 

and  /i,/2,  /a  are  2D  harmonic  functions.  Furthermore,  it  may  be  noted  that  the  first  series  has 
complex  eigenvalues  and  eigenfunctions  while  the  second  has  only  real  eigenvalues  and 
eigenfunctions.  For  an  explicit  definition  of  all  functions  see  Penado  and  Folias  (1989). 

4.  LOADING  TRANSVERSE  TO  THE  FIBERS 

( i)  Interior  Stress  Field: 

For  a  uniform  transverse  loading  along  the  fiber  direction  and  under  the  assumption  that 
perfect  bonding  prevails  at  the  fiber  /  matrix  interface,  the  3D  stress  field  (at  the  interface)  and 
along  the  fiber  length  is  found  to  be  constant  (  see  Fig.  2 )  all  along  the  interior  and  that  as  one 
approaches  the  free  surface  a  boundary  layer  is  noted  to  prevail  where  the  stress  field  increases 
rather  rapidly.  This  rapid  change  suggests,  therefore,  the  presence  of  a  possible  stress  singularity. 
Moreover,  the  width  of  this  3D  boundary  layer  is,  approximately,  two  fiber  diameters  from  the 
free  edge.  The  reader  may  also  note  that  a  second,  3D,  effect  is  that  the  amplitude  of  the  stresses 
at  the  center  of  the  fiber  length  is,  in  general,  a  function  of  the  ratio  of  fiber  diameter  /  fiber 
length.  If  that  ratio,  however,  happens  to  be  less  than  or  equal  to  1/10,  then  all  along  the  interior  a 
'pseudo  plane  strain'  condition  prevails.  Figs.  2a,  2b  depict  the  profile  of  the  interface  stresses  at 
0=0  as  a  function  of  z  /  h.  The  numerical  results  are  specialized  for  the  material  system  :  SCS-6  / 
fibers,  Beta21  /  matrix.  Figs.  3,  and  4,  show  typical  interface  stress  profiles  of  the  matrix  and  fiber 
on  the  plane  z=0,  and  as  functions  of  the  angle  0.  The  reader  may  notice  that  the  max.  of  the 
stresses  occurs  at  the  location  0=0. 

(ii)  Edge  stress  Field: 

As  it  was  previously  noted,  in  the  neighborhood  where  one  approaches  a  free  surface  e.g. 
the  edge  of  the  plate  or  in  the  vicinity  of  crack  bridging  ( see  Fig.  5  below),  there  may  very  well 
be  present  a  stress  singularity.  Utilizing  a  local,  3D,  asymptotic  analysis  one  can  substantiate  the 
presence  of  a  weak  stress  singularity.  Complete  details  of  this  analysis  may  be  found  in  the  work 
of  Folias  (  1989  ).  Without  going  into  the  mathematical  details,  a  summary  of  the  results,  at  room 
temperature  and  for  the  composite  material  system  discussed,  is  given  below: 
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NORMALIZED  STRESSES  AT  THE  INTERFACE 

MATERIAL  SYSTEM  :  SCS-6  /  FIBER,  BETA21  /  MATRIX 

LOADING  :  UNIFORM  TENSION  TRANSVERSE  TO  THE  FIBER 
FIBER  DIAMETER  /  FIBER  LENGTH  =  1/10 


z/h 


Fig.  2a.  Interface  matrix  stresses  as  a  function  of  z/h. 
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FREE  EDGE 


NORMALIZED  OCTAHEDRAL  SHEAR  STRESS  AT  THE  INTERFACE 


MATERIAL  SYSTEM  :  SCS-6  /  FIBER,  BETA21  /  MATRIX 

LOADING  :  UNIFORM  TENSION  TRANSVERSE  TO  THE  FIBER 
FIBER  DIAMETER  /  FIBER  LENGTH  =  1/10 


Fig.  2b.  Interface  matrix  octahedral  shear  stress  as  a  function  of  z/h. 
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NORMALIZED  STRESSES  AT  THE  INTERFACE 

MATERIAL  SYSTEM  :  SCS-6  /  FIBER,  BETA21  /  MATRIX 

LOADING  :  UNIFORM  TENSION  TRANSVERSE  TO  THE  FIBER 
FIBER  DIAMETER  /  FIBER  LENGTH  =  1/10 


CORNER  ANALYSIS  AT  THE  EDGE 


Possible  composite  failure  mode 


FIBER  MEETING  A  FREE  EDGE 

Local  3D  Stress  Field  :  ( Folias  IJF  1989 ) 

a,;  =  ())) 

where  for  a  Titanioum  matrix  and  SiC  fibers 
a  =  0. 1 10,  at  room  temperature 
a  =  0.190,  at  900®  C 

( i)  location  :  (j)  =  ^,for  Gf/Gm  =  3.608 

o^”^  =  -11.219p^B^ 

G^ee^  =  -4.823p-“5W 

(ii)  location  :  <}>  =  0. 

0^rr^=-9.97p-“BW 
G^e9^  =  -3.391p-“5<'">. 

Fig.  5.  Geometrical  configuration  and  basic  results. 
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OCTAHEDRAL  SHEAR  STRESS  AT  THE  CORNER  OF  A  FREE  EDGE; 


REMARK  :  Crack  will  initiate  at  the  fiber  edge  and  at  the  interface  and  will  then 
follow  the  direction  of  the  maxnretahedral  shear  stress. 


Fig.  5b.  The  octahedral  shear  stress  at  the  comer. 


The  following  observations  are  worthy  of  note.  First,  as  the  ratio  of  the  shear  moduli 
increases,  so  is  the  stress  singularity.  This  is  compatible  with  our  physical  expectations  and  within 
the  assumptions  of  our  theory.  Second,  all  things  being  equal,  at  the  edge  the  controlling  stress  for 
failure  is  the  radial  stress  particularly  at  the  location  (f)=0  and  0=0. 

Room  Temp.  900^  C 

at  (t)  =  0  ;  a{?’Vo^e^=  2.94  =  2.94 

at<l>  =  7t/2:  ar/a^^  =  2.33  ar/a^'=2.18 

Similarly,  in  the  vicinity  of  the  edge,  the  octahedral  shear  stress  attains  a  maximum  at  an 
angle  <t)=40  degrUs  (see  Fig.  6).  It  is  interesting  to  note  that  in  this  neighborhood,  the  ratio  of 

Xocfmax  Iroomfemp/Xocfrnax  I900°C  =  1  /2.09 

This  suggests,  therefore,  that  as  AT  increases  the  application  of  a  transverse  loading  will  cause  the 
matrix  to  undergo  substantial  more  plastic  deformation  in  this  region. 

Computing  next  the  displacement  at  the  free  surface  z  =  h,  we  notice  that  its  magnitude 
increases  at  elevated  temperatures. 

wl^,0=O,/?oomremp/W'4=O.e=O.9CX3OC  =  1/1  -6. 

Physically,  this  suggests  that  a  mode  I,  H  and  IH  crack  failure  will  initiate  at  the  edge  and 
at  the  interface  due  to  an  applied  load  transverse  to  the  fibers. 
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5.  RESIDUAL  STRESSES  DUE  TO  A  THERMAL  LOAD. 


Residual  stresses  due  to  curing  and  thermal  stresses  due  to  differences  between  the 
thermal  expansion  coefficients  of  the  matrix  and  fiber  may  have  a  major ’effect  on  the 
microstresses  within  a  composite  material  system  and  must  be  added  to  the  stresses  induced  by 
the  external  mechanical  loads.  Such  microstresses  are  often  sufficient  to  produce  microcracking 
even  in  the  absence  of  external  loads,  example  during  the  cooling  process. 

In  this  investigation,  a  systematic,  3D,  micromechanics  approach  is  used  in  which  the , 
fibers  of  a  composite  material  system  are  modeled  as  cylindrical  inclusions  which  are  embedded 
into  a  matrix  plate.  The  analytical  model  is  then  used  to  predict,  the  residual  stresses!  due  to  a 
thermal  expansion  mismatch,  eg.  during  the  cooling  process.  Moreover,  the  model  providesoa 
better  understanding  of  how  the  residual  stresses  are  developed  and  how  they  can  be  conhEplfed 
particularly  in  relation  to  ceramics  where  there  is  no  ductility  to  accommodate  plastic 
deformation. 

The  analysis  reveals  the  dependence  of  the  residual  stress  field  on  the  fiber  volume  fraction 
ratio,  identifies  the  critical  locations  where  a  crack  is  most  likely  to  initiate,  recovers  the  interface 
shear  stress  profile  and  provides  important  information  to  the  material  designers  for  the 
pre-selection  of  fiber  and  matrix  materials  in  order  to  alleviate  some  of  the  residual  stresses. 

Without  going  into  the  mathematical  details,  we  consider  a  composite  material  system 
consisting  of  SIC-6  fibers  which  are  embedded  into  a  beta21  matrix  plate  and  the  entire  system  is 
then  exposed  to  an  environment  of  a  uniform  cooling  temperature  AT.  While  it  is  true  that  the 
material  constants  do  change  as  a  function  of  the  temperature,  the  thermal  coefficients  appear  in 
the  solution  as  a  ratio  and  interestingly  enough  this  ratio  changewary  little.  On  the  other  hand,  the 
ratio  of  the  shear  moduli  changes  considerably  as  the  temperature  varies.  Thus,  the  results- are 
very  much  dependent  on  the  material  properties  which  one  uses.  Thus,  if  one  bases  the  analysi?  on 
the  shear  moduli  ratio  at  room  temperature,  the  following  stress  profiles  are  recovered  at  :he  fiber 
/  matrix  interface.  Fig.  7  depicts  the  radial  matrix  stress  on  the  plane  z=0  and  as  a  function  of  the 
angle  0  .  It  is  noted  that  the  radial  stress  is  compressive.  Similarly,  the  tangentiaT  stress  is  tensilein 
nature  and  its  maximum  occurs  at  the  location  of  0  =  0.  In  general,  the  location“of  this  maximum;, 
is  a  function  of  the  material  properties  and  particularly  the  shear  moduli  ratio.  Moreover,  in  the 
above  analysis  perfect  bonding  was  assumed  to  prevail  at  the  fiber  /  matrix  interface.  If,  however, 
we  relax  the  conditions  at  the  interface  and  allow  slippage  then  the  maximum  occurs  elsewhere. 
More  specifically  in  this  case  it  occurs  at  0  =  45. 

Examining  next  the  possibility  of  matrix  cracking,  it  becomes  evident  from  the.abpve  that 
no  cracking  will  occur  in  the  matrix  for  a  AT  =  900°C.  This  matrix  material  is  too  strpng  for 
preexisting  microcracks  to  grow.  Examination  of  the  azz  stress  also  shows  that  no  cracks  will 
develop  in  that  direction  either.  These  results  are  in  line  with  the  obtained  in  house  fesults  based 
on  finite  elements  ( J.  Kroupa,  1994). 
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material  system  :  SCS-6  /  FIBER,  BETA21  /  MATRIX 
loading  :  UNIFORM  THERMAL  LOADING 


7.  Interface  radial  matrix  stress  on 


the  plane  z=0  and  as  a  function  of 


INTERVACi:  TANGENTIAL  S  I'KESS 


MATERIAL  SYSTEM  ;  SCS-<»  /  FIBER,  BETA21  /  MATRIX 
LOADING  :  UNIFORM  THERMAL  LOADING 


Fig.  8.  Interface  aee  matrix  stress  on  the  plane  z=0  and  as  a  function  of  9. 
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Finally,  it  should  be  noted  that  the  effect  of  the  shear  moduli  ratio  on  the  interface  stresses 
is  substantial.  This  can  be  seen  by  the  following  comparison  of  the  tangential  interface  stress  when 
using  three  different  shear  modulii  ratios  which  reflect  different  temperature  levels; 


a,3  /(a„G„  AT)  /(a^G^  AT)  /(a^G^  AT) 

_V  Room  Temp.  Mid  Temp.  HiehTemp. 

0  39  1.36  1.59  1.64 

Although  it  would  be  desirable  to  have  a  program  in  which  the  material  properties  do  vary  with 
temperature,  one  can  compensate  by  taking  the  results  corresponding  to  the  high  shear  modulii 
ratio.  The  thermal  expension  coefficients  on  the  other  hand  appears  as  a  ratio  which  ratio  does  not 
vary  appreciably  to  make  any  significant  differences. 

The  variation  of  the  normalized  tangential  interface  stress  as  a  function  of  the  fiber  volume 
fraction,  for  this  material,  is  almost  linear  and  may  be  approximated  with  the  equation 

-?a-=  =  0.92+1.02K/+0.28  fj.  (13) 


6.  CONCLUSIONS. 

In  view  of  the  above,  the  matrix  will  not  exhibit  any  cracking  as  a  result  of  the  residual 
stresses  which  are  developed  during  the  cooling  process  (  AT=  900®C).  Moreover,  the  residual 
stresses  predicted  are  in  very  good  agreement  with  those  obtained  in  house.  (  WL/MLLN  ) 
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1  Introduction 


In  design,  anticipation  of  material  failure  is  vital.  Over  the  last  decades,  a 
theory  first  proposed  by  Griffith  has  become  one  of  the  most  used  tools  in 
prediction  of  fracture  of  simple  materials.  In  the  Griffith  theory,  a  material 
is  assumed  to  have  microscopic  cracks  with  high  stresses  found  at  their  tips; 
these  cracks  cause  the  material  to  fail  at  a  much  lower  level  than  molecu¬ 
lar  binding  forces  predict.  To  accurately  predict  failure  of  these  materials  re¬ 
quires  knowledge  of  the  stress  field  near  the  tips. 

Composite  materials  also  have  cracks,  but  models  are  much  more  involved 
because  of  material  interactions.  As  before,  accurate  determination  of  the 
stress  field,  especially  near  high  stress  regions,  is  vital  for  failure  prediction. 

For  all  but  the  simplest  geometries,  closed  form  solutions  are  almost  impos¬ 
sible  to  obtain,  so  numerical  methods  have  become  very  popular.  For  ma¬ 
terial  science,  by  far  the  most  popular  methods  are  the  finite  element  (FE) 
methods.  Some  of  their  desirable  properties  are  relative  conceptual  simplic¬ 
ity,  straightforward  (but  tedious)  implementation,  a  large  available  code  base 
for  solving  problems,  and  a  tendency  to  require  only  modest  computing  re¬ 
sources.  Since  finite  elements  are  based  on  polynomials,  the  solutions  to  prob¬ 
lems  for  which  finite  elements  are  used  must  be  expressible  as  a  sum  of  poly¬ 
nomials,  or  the  FE  approximation  will  be  very  poor.  Since  functions  wifii 
imbounded  derivatives  (i.e.  stresses  near  cracks  and  material  interfaces)  can¬ 
not  be  approximated  well  with  polynomials,  it  is  common  to  use  special  el¬ 
ements  for  these  functions. 

For  example,  in  2  dimensional  problems  and  some  special  3  dimensional  cases, 
it  has  been  shown  that  the  stresses  near  the  crack  tips  are  proportional  to 
1  /  -y/r,  and  that  the  constant  of  proportionality,  Kc,  depends  only  on  the  ge¬ 
ometry  of  the  material.  This  is  enough  information  to  complement  the  (poly¬ 
nomial  based)  FE  methods  with  singular  elements  (which  behave  like  l/\/r 
in  the  appropriate  regions),  thus  reducing  the  problem  to  one  for  which  FE 
are  well  suited,  and  get  good  numerical  answers. 

For  other  3  dimensional  crack  problems  however,  asymptotic  expansions  have 
shown  the  stress  fields  to  be  proportional  to  r““,  0  <  a  <  1/2,  a  depending 
on  the  geometry  of  the  material.  Thus,  for  3D  problems,  one  does  not  in  gen¬ 
eral  know  the  behavior  of  the  singularity  a  priori,  and  the  problem  cannot 
be  reduced  to  one  which  FE  can  handle  well,  resulting  in  low  accuracy  near 
singularities.  Notice  the  difficulty  here:  the  stress  distribution  in  the  vicinity 
of  a  crack  also  depends  strongly  on  material  geometry,  making  it  impossible 
to  separate  behavior  near  cracks  from  the  rest  of  the  material —  but  this  sep¬ 
aration  is  how  FE  methods  handle  these  problems. 

This  inherent  weakness  of  most  numerical  mefirods  —  the  inability  to  han¬ 
dle  singularities  without  "assistance"  —  is  not  shared  by  the  group  of  sinc- 
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function  based  methods.  Further,  sine  methods  enjoy  an  exponential  con¬ 
vergence  rate,  enabling  one  to  get  many  digits  of  accuracy  with  reasonable 
work,  if  desired. 

For  crack  and  related  problems,  this  means  only  the  location,  but  not  the  type, 
of  the  singularity  is  needed,  and  the  solution  can  be  accurately  computed. 

The  weakness  of  using  these  methods  lie  in  the  fact  that  their  use  for  differ¬ 
ential  and  integral  equations  is  recent,  so  no  large  code  base  yet  exists;  their 
use  has  been  largely  in  ID  problems;  and  they  are  not  as  intuitive  as  FE  or 
finite  difference  methods,  hence  overlooked  by  most  people. 

The  purpose  of  this  work  over  the  last  years  has  thus  been  the  further  devel¬ 
opment  of  sine  methods  for  systematic  solution  of  problems  in  mechanics 
that  possess  singularities.  The  development  of  the  method  has  been  guided 
by  a  representative  problem,  described  in  detail  below. 

It  should  be  emphasized  that  the  following  problem  is  only  one  example 
of  an  entire  class  of  problems,  serving  here  to  illustrate  the  effectiveness 
and  flexibility  of  the  sine  methods  for  that  class  of  problems. 

It  should  also  be  noted  that  this  method  is  readily  extended  to  handle  fully 
nonlinear  material  behavior. 

The  details  of  the  method  and  development  to  date  comprise  the  bulk  of  the 
remainder  of  this  report. 


2  Problem  and  Solution  Approach 

The  problem  here  is  to  find  the  stress  field  uniformly  to  a  desired  accuracy^ 
in  the  piece  of  composite  material  shown  in  figure  1.  To  this  end,  we  solve 
the  full  isotropic  Navier's  equations 

V"u  +  j-^V(V  .  u)  =  0  (1) 

with  appropriate  displacement  and  stress  boundary  conditions. 

This  problem  is  broken  down  as  follows. 

1.  Using  this  problem's  radial  S5anmetry  allows  immediate  reduction  to 
a  sequence  of  2D  problems  on  domains  as  shown  in  figure  2  —  section 
2.1. 

2.  At  points  Ai,...,  A4,  singularities  are  to  be  expected.  Typically,  these 
singularities'  behavior  depends  on  approach  direction,  so  each  of  these 
problems  is  further  divided  into  triangles,  and  these  triangles  are  mapped 
to  coupled  rectangles  as  shown  in  figures  3  and  4  —  section  2.2. 

^Uniform  accuracy  of  3  digits  in  both  displacements  and  stresses  is  easily  achieved. 
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Figure  1:  The  domain  for  the  full  3D  problem.  The  debonded  region  is 
treated  as  a  very  thin  crack. 


3.  The  sine  method  is  applied  to  the  resulting  collection  of  rectangles,  and 
the  solution  is  thus  obtained  —  section  2.3. 

The  details  of  steps  (1)  and  (2)  are  currently  being  worked  on,  and  are  pro¬ 
gressing  rapidly.  We  therefore  only  present  an  outline  of  the  required  steps 
in  sections  2.1  and  2.2.  The  sine  method  and  its  application  to  (3)  are  found 
in  section  2.3. 


2.1  Forming  of  the  2D  sequences 


Although  the  sine  methods  can  be  used  directly  for  3D  problems,  incremen¬ 
tal  development  is  best  done  for  2D  problems  first.  For  the  given  problem, 
we  take  advantage  of  the  geometry  and  use  the  substitutions 


Ur{T,e,z) 
we(r,  d,  z) 
u^{r,e,z) 


oo 

Urfi{r,  2)  +  13  “r,n(r,  z)  COS 
n=l 


00 


13  z)  sin 

n=l 


00 


Wz.o(r,  ^)  +  13 

n=l 


(2) 

(3) 

(4) 
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Figure  2:  Reduced  3D  problem.  Crack  and  points  with  singularities  are  in¬ 
dicated. 


in  the  polar  form  of  equation  1  to  reduce  the  original  equations  to  the  men¬ 
tioned  sequence  of  2D  problems.  The  equations  governing  figure  2  are  thus 


52 

dz^ 


Uz,0 


I 


-f- 


52 


\dzdr 


«r,0  =  0 


(6) 


for  the  non-0  terms  of  equations  2-4  and 


-2{v 


-1) 


[dr^  “^’7  "^2  r2 


(-1  +  21/ 


r 


2  Ur,n  1  {-l+2u)  Ur, nU^ 
7.2  4  J.2 


+ 


(7) 
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Figure  3:  Triangulation  of  reduced  problem. 


ip  ^2 

1  (u-l)ue,nn^  1  i-3  +  4u)ur,nn 

2  2 

_  (-1  +  2t/)  ^  1  ^r,n)  ^ 

T  2  7*^  7*^ 

1  _  n 

r3  2  r2  ” 


(8) 


-2(i^-l) 


(-1  +  21/ ) 


I  ^  (-l  +  2i/)«;,,„n2 
^'"’7  "^4  r2 


(  g"  \  ,  1  ^  _  (-1  +  21/)  (|:t/.,n)  ^ 

\dzdr^'^''^)  2  r  r  r 


=  0 

for  the  9  terms  -  notice  the  n  dependence  here. 


(9) 


So  far,  we  thus  have  one  2  imknown  system,  and  one  3  unknown  system, 
both  independent  of  6.  Of  course,  each  set  occurs  twice — once  for  the  matrix 
and  once  for  the  fiber,  and  appropriate  coupling  is  done  across  the  botmd- 
aries  by  matching  stresses  and  displacements  above  and  below  the  crack, 
and  requiring  zero  stresses  at  the  crack. 

The  boimdary  conditions  are  transformed  similarly. 
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Figure  4:  Rectangular  solution  regions  from  triangulation  of  reduced  prob¬ 
lem.  Notice  that  the  four  domains  shown  share  the  boundary  that  came  from 
the  point  Ai 


2.2  Mapping  of  the  Triangulation 

After  the  reduction  to  a  sequence,  there  may  be  singularities  at  the  points 
Ai, Ai  (figure  3).  Most  likely,  these  will  be  moving  singularities,  e.g.  us¬ 
ing  a  second  polar  coordinate  system  with  origin  at  point  Ai,  the  displace¬ 
ments  may  have  a  form  similar  to  /  =  r“  sin(0),  0  <  a  <  1/2.  In  the  coor¬ 
dinate  system  used  in  figure  2,  /  would  behave  like  z{r^  -f-  22^0/2-1/2 
Ai,  and  its  partials  would  have  the  dominant  term  rzf{T^  +  22)2/2-0/2  jj^g 
form  of  singularity  is  easily  changed  to  a  more  suitable  form  via  the  map¬ 
pings  r  =  z  =  Geometrically,  this  means  splitting  rectangular  do¬ 
mains  into  triangular  pieces,  as  shown  in  figure  3.  Algebraically,  tiiis  map¬ 
ping  turns  dominant  terms  of  derivatives  into  the  form  which  no 

longer  depends  on  approach  direction,  enabling  direct  solution  using  sine 
methods. 


2.3  Sine  Method  Description 

It  is  the  purpose  of  this  chapter  to  give  necessary  background  for  the  under¬ 
standing  of  sine  methods  as  used  for  the  present  work.  We  begin  with  some 
known  one-dimensional  results  in  section  2.3.1;  the  extension  to  two  dimen¬ 
sions  is  shown  in  section  2.3.2.  The  extension  to  higher  dimensions  is  similar, 
and  will  not  be  shown  here. 

Section  2.3.2.5  demonstrates  a  2  dimensional,  multiple  unknown,  multiple- 
domain  problem  of  the  type  occurring  after  the  triangulation  done  above. 
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2.3.1  Sine  methods  for  ID  problems 

In  this  section,  we  summarize  some  available  results  for  sine  interpolation 
and  sine  collocation.  The  main  reference  for  technical  details  of  this  section  is 
the  book  Numerical  Methods  based  on  Sine  and  Analytic  Functions,  Springer, 
1993;  most  results  are  proven  there. 


2.3.1.1  Interpolation  and  simple  collocation  First,  some  definitions. 
The  sine  (or  Whittaker  Cardinal)  function  is  defined  by 

sin(7r  x) 


sinc(a:)  = 


'KX 


(10) 


Define  the  domain  Dd  by 

Dd  =  {we  C\^{w)  <  d}  (11) 

Let  a  >  0,  and  let  La  (Dd)  denote  the  family  of  functions  /  with  the  properties 

•  /  is  analytic  in 

•  for  some  c  >  0  and  all  2;  G  Dd, 

^OCZ 

l/(^)l  <  (12) 


Now,  taking  h  = 


^aN) 


'(l  +  |e*|)2« 

,  we  have  the  interpolation  result 


nx)=  E  /(*fe)sinc(^-^)+E(W) 

E{N)  =  ciVNe-^^^^  (13) 

for  a  positive  ci  depending  only  on  /,  d  and  a.  Notice  that  this  result  is  for 
the  real  line,  and  the  function  must  decay  at  ±00. 

By  first  remapping  functions  approaching  a  nonzero  limit,  this  can  be  en¬ 
hanced  to  handle  non-zero  values  at  ±00: 

f(x)=  E  c,smc(^^)+c„+iS<„(a:)  +  c_„-,S-„(i)  +  £:(itf)  (14) 

k=-M  \  "  / 


Ck  =  f{kh),  k  =  -N..N 

(15) 

II 

1 

1 

(16) 

CN+l  =  /(OO) 

(17) 

5-00  ((c)  =  - 

^  1  -1-  e®® 

(18) 

^OeX 

Soo(x)  = 

^  '  1  -1-  e®® 

(19) 
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Notice  that  the  summation  runs  from  -M  to  A^;  the  error  bound  is  of  the 
same  form  as  above. 

Defining  masm  =  2A^+ 1,  dus  means  once  n  digits  of  accuracy  are  obtained, 
one  can  roughly  get  I  An  digits  by  doubling  m. 

To  interpolate  a  function  /  defined  on  [a,  6]  c  R  we  first  make  the  following 
definitions: 

1.  for  [a,  6]  €  D,  ^  is  a  conformal  map  with  4>’-  D  -¥  Dd 

2.  ^  =  <f>-^ 

3.  p  = 

4.  =  sine 


Let  a  >  0,  and  let  La  (D)  denote  the  family  of  functions  F  with  the  properties 


•  F  is  analytic  in  D 

•  for  some  c>  0  and  aHz  e  D, 


W|  <  c 


p(^)" 

(1  +  \p(z)\)‘^ 


(20) 


(7rd  \  ^  ^ 

j  ,  we  have  the  interpolation  result 
f{^)=  Z!  Ck'ik'\^) E{N) 

kz=^N 

E{N)  =  Cy/Ne-'^^^  (21) 

for  a  positive  ci  depending  only  on  /,  d  and  a.  Notice  that  equation  20  re¬ 
quires  /  to  vanish  at  the  endpoints  of  the  interval. 

As  before,  this  series  can  be  enhanced  to  handle  nonhomogeneous  endpoint 
values  with  a  simple  addition: 


F'ix)=  £  Cfc7W(a:)-t-c^+i56(a;)-H.c_M-i5a(:r)-hF(M,iV)  (22) 

Ck  =  {F-Sa-  Sb){i>{kh)),  k  =  -N..N  (23) 

cn+1  =  F{b)  (24) 

c-M-i  =  F{a)  (25) 
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and  Sa  and  Sb  are  cubic  splines  with  value  1  at  the  left  and  right  endpoints, 
respectively,  and  zero  derivatives  at  the  endpoints. 

As  written,  the  expressions  for  the  c*  are  no  longer  simple  function  evalua¬ 
tions.  This  extra  work  can  be  shifted  to  the  expansion  of  F  by  defining  the 
discrete-orthogonal  terms 


Sa{x)  =  Sa{x)  -  dklk{x)  (26) 

k=-M 

dk  =  Samh))  (27) 

Sb{x)  =  Sb{x)-  Yl  dkjkix)  (28) 

k=-M 

dk  =  Sbmkh))  (29) 

With  these  definitions,  we  have  the  expansion 

■^(^)  =  H  CA:7f^(a:) +C65ft(a;) +  Co5a(a;)  (30) 

k=-N 


with  Ck  =  F{xk),  Ca  =  F{a),  ct  =  F{h),  and  Xk  =  il){kh) 


2.3.1.2  Extensions  for  mixed  boundary  value  problems  For  the  solution 
of  mixed  boimdary  value  problems,  a  finite  nonzero  approximation  of  the 
derivative  at  endpoints  is  needed.  Since  the  derivatives  of  the  '^k  are  im- 
boimded  at  the  endpoint,  a  nullifier  g  is  introduced  to  make  the  derivatives 
of  the  series  terms  also  vanish  at  the  endpoints.  Then,  as  before,  adding  extra 
terms  with  the  right  properties  gives  a  useful  basis. 

Let  Ta{x)  and  Tb{x)  be  cubic  splines  with  derivatives  of  one  at  the  left  and 
right  endpoint,  respectively,  and  other  values  and  derivatives  zero  at  the  end¬ 
points;  then  the  following  series  can  be  used  to  approximate  /  on  [a,  6]  when 
/  is  specified  via  a  mixed  boimdary  value  problem. 


N 


Y  ^k7k(x)  +  CbSb(x)  +  CaSa(x)  +  Ci/fb(x)  + 
k=-M 

(31) 

Ca>t(x)  +  E(N,M) 

(32) 

.  ( <j){x)  -  kh\  , 

7fc  =  smcl  ^  \g{x) 

(33) 

Xk  —  i’ikh) 

(34) 

II 

(35) 

Ck  =  f{xk)/9{xk) 

(36) 

9 


Cb 

Ca 

cv 

Oq/ 

Sa{x) 


dk 


dk 


faix) 


dk 


fb(x) 


dk 


=  m 

=  /(a) 

=  m 

=  f'ia) 

=  Sa{x)-  dkjkix) 

k=-M 

SgjXk) 

gi^k) 

=  Sb{x)-  dknfkix) 

k=-M 

_  Sbjxk) 

9{xk) 

N 

=  Ta{x)  -  dkjkix) 

k=-M 

_  Ta{Xk) 

9{xk) 

N 

=  Tb{x)  -  dkjkix) 

k=-M 

Tbjxk) 

9{xk) 


(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 


All  following  extensions  are  based  upon  this  series  representation,  or  sub¬ 
sets  of  it. 


2.3.2  Sine  methods  extended  for  2D  problems 

The  straightforward  way  to  extend  the  sine  series  to  higher  dimensions  is 
via  tensor  product.  In  this  section,  two  derivations  of  the  series  expansion 
are  given  and  the  matrix  structure  which  arises  for  collocation  is  described. 

In  section  2.3.2.1,  we  first  show  the  derivation  under  the  assumption  that  the 
unknown  can  be  fully  represented  by  the  series  (31),  and  using  a  tensor  prod¬ 
uct  to  get  the  two  dimensional  extension. 

Next,  we  treat  the  splines  in  (31)  as  a  change  of  unknown,  so  as  to  produce  a 

N 

homogeneous  equation  to  be  satisfied  by  the  simple  Stun  ^  Ckjk{x),  and 

k=-M 

extend  this  line  of  thinking  to  two  dimension  -  section  2.3.2.2. 

Extension  of  the  series  to  multiple  domains  is  considered  in  section  2.3.2.3 

Splitting  of  the  two  dimensional  series  into  logical  units  for  collocation  and 
matrix  setup  is  the  subject  of  section  2.3.2.4. 
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2.3.2.1  Tensor  product  derivation  Let  u  {x)  be  represented  by  the  full  sum 
used  in  equation  (31),  written  as  single  sum.  Thus 

N+2 

=  'll  Cfc7fe(a;)  =  Ckjk{x)  (49) 

k=-M-2 

where  the  last  equation  is  written  in  summation  notation  (repeat  subscripts 
are  summed  over).  Then  we  have 

u{x,y)  =  [dkiiy)]lki{x) 

=  [4ifc27fe2(2/)]7iki(a:) 

Ni+2  JV2+2 

=  £  £  dk^k2ikx{x)ikAy)  (50) 

ki=-Mi-2  k2=-M2-2 

This  representation  is  valid  on  any  rectangular  region;  regions  with  other 
shapes  can  easily  be  mapped  onto  a  rectangle.  Also,  the  "rectangle"  can  be 
imbounded  on  one  or  more  sides;  only  the  choice  of  conformal  map  (below) 
changes.  Thus,  half-  or  fullspace  problems  can  be  solved  easily. 

Assume  for  simplicity  Ni  =  N2,  Mi  =  M2  .  Define  masm  =  N  +  M  +  l. 
Then  the  series  (50)  has 

(m  -1-  4)^  =  -I-  8  m  -h  16  (51) 

terms.  Notice  that  this  is  the  most  general  form  possible,  used  for  problems 
with  mixed  conditions  on  all  boimdaries.  Since  not  every  problem  has  mixed 
conditions  on  all  boimdaries,  the  nullifier  must  be  selected  in  conjunction 
with  the  splines  for  each  direction  and  boundary  separately  before  forming 
the  tensor  product. 

To  this  end,  it  is  more  natural  to  think  of  the  splines  as  a  remapping  of  the 
unknown  -  the  subject  of  the  next  section. 

For  collocation  points,  starting  from  {xk\xk  =  U  {xa,Xb} 

for  the  one-dimensional  case,  we  get  the  collocation  points  as  a  tensor  prod¬ 


uct  also: 

{Xkj\Xkj  =  {'ll’{khi),t{){jh2))}k=-Mi-l..Ni+lJ=-M2-l..N2+l  U  (52) 

{X(ij\Xaj  —  (^^aj  ^(j/l2))}'j=— Afj— I..N2+I  ^  (5^) 

{Xbj\Xbj  =  {Xb,‘tp{jh2))}j=-M2-l..N2+l  U  (54) 

{ajfeakjka  =  ii’{khi),ya)}k=-Mi-l..Ni+l  U  (55) 

{xkb\xkb  =  {‘il^{khi),yb)}k=-Mi-i..Ni+i  (56) 


Notice  that  we  have  (m  -t-  4)^  =  -t-  8  m  -I- 16  points,  as  expected. 

2.3.2.2  Tensor  products  revisited  Recall  the  simple  series  used  for  La  func¬ 
tions  in  ID  (equation  (21)).  By  using  a  nullifier,  this  series  has  zero  value  and 
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zero  derivative  at  the  boundaries.  By  providing  splines  for  the  nonhomoge- 
neous  boundary  conditions  (equation  (22)),  the  problem  Lu  =  f,Bu  =  g  is 
effectively  remapped  to  Lv  -  f,Bv  =  0  and  this  problem  can  be  solved  with 
the  sine-only  series  (21). 

The  same  approach  can  be  taken  in  2  (and  higher)  dimensions.  Starting  with 
the  simple  sine  series  and  forming  the  tensor  product,  we  obtain  the  repre¬ 
sentation 


Ni  N2 

y)  =  £  £  ^i)  °  {^)  {SU,  ^2)  o  h92)  {y)  (57) 

jzz^Mi  A:=— M2 

valid  for  u  €  La,  u'{a)  =  u'{b)  =  0.  It  thus  remains  to  remap  the  problem 
Lu  =  f,  Bu  =  g  to  Lv  =  f^Bv  =  0.^  Taking  a  hint  from  the  previous  section, 
on  each  boimdary  of  the  rectangle,  we  can  use  a  series  of  the  form 

j=-Mi 

for  u's  value  and  tangential  derivative,  and  a  series  of  the  form 

Z!  [('5'(^.  h)  O  ^g)  (xi)  -h  5a(rri)  -H  56(a;i)l  f{x2)  (59) 

j=-Mi 

for  1^.  The  direction  of  the  boundary  is  xi,  the  normal  direction  X2.  This 
gives  a  total  of  -h  8m  H-  24  terms  per  unknown.  The  causes  for  this  larger 
number  are  redimdancies  at  the  comers  of  the  domain,  in  the  spline-only 
terms.  First,  u  only  has  one  value  at  each  comer,  reducing  the  total  number 
of  terms  by  4.  Then,  the  x  and  y  partials  at  each  comer  are  unique  also,  fur¬ 
ther  reducing  the  number  of  terms  by  8.  Thus,  we  are  left  with  m^  +  8m  + 12 
terms  per  unknown.  The  four  terms  missing  here  but  present  in  (50)  are  the 
second  order  mixed  derivative  terms  representing  at  the  comers  of  the 
domain. 

The  results  in  sections  2.3.2.1  and  2.3.2.2  thus  differ  at  the  comers  of  the  do¬ 
main;  numerically,  the  presence  of  terms  representing  the  mixed  second  par¬ 
tial  is  redundant  and  these  terms  are  not  used. 

Next,  we  derive  at  the  collocation  points  by  extending  the  ID  problem  in  a 
way  analogous  to  the  series  derivation,  rather  than  by  tensor  product.  First, 
recall  that  for  ID,  tiuree  point  regions  can  be  distinguished.  Always  present 
are  interior  collocation  points,  given  by  x*  =  'ip{kh),  k  =  -M..N  since  these 
correspond  to  the  sine  part  of  the  series.  When  the  splines  representing  the 

^Note:  for  single  imknown  problems,  an  explicit  remapping  can  usually  be  found,  so 
that  the  resulting  problem  has  only  the  coefficients  of  (57)  as  unknowns.  However,  with 
multiple  unknowns,  only  the  form  of  the  remapping  is  known,  and  unknown  coefficients  of 
this  form  also  go  into  the  matrix  system.  Thus,  this  approach  is  somewhat  implicit. 


[{S{k,  h)  O  01^1)  (xi)  -b  5a(a;i)  -b  56(xi)  -b  fa(xi)  -b  ffe(xi)]  S(X2)  (58) 
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Figure  5:  The  point  layout  for  2D  collocation,  for  N  =  M  =  1.  The  circles  are 
the  boundary  points,  while  the  crosses  represent  the  always  present  part  of 
the  interior  points  obtained  from  the  ID  tensor  product.  The  points  marked 
with  triangles  are  used  only  when  the  derivative  spline  terms  corresponding 
to  them  are  present.  Notice  that  the  points  are  not  really  evenly  spaced. 
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value  of  the  unknown  at  the  endpoints  of  the  interval  are  present  (5x„  etc.), 
the  collocation  point  set  is  expanded  to  include  the  boimdary  points.  Lastly, 
when  the  splines  representing  the  derivatives  at  endpoints  (Tx^  etc.)  are  in¬ 
cluded,  the  collocation  points  are  further  extended  with  extra  points  in  the 
interior:  Xk  =  'ip{kh),  k  =  {—M  —  1,N  +  1} 

Selecting  the  collocation  points  in  2D  analogously  we  have  a  point  grid  as 
shown  in  figure  5. 

Notice  that  here  we  have  no  points  at  the  comers,  unlike  equations  (52);  this 
is  consistent  with  the  series  selection  in  equations  (57)  through  (59),  since 
there  we  also  have  4  terms  fewer  than  equation  (50). 

2.3.2.3  Multiple  Domain  problems  in  2D  Handling  of  multiple  domain 
problems  is  straightforward.  At  the  connecting  boundaries,  the  equations 
involve  unknowns  from  both  domains  and  this  simply  reflects  in  the  collo¬ 
cation  matrix. 

Stated  another  way,  the  continuous  properties  of  the  linear  problem  are  re¬ 
flected  in  the  discrete  approximation  via  the  linear  system.  Therefore,  one 
has  to  only  consider  the  matrix  implications  of  domain  coupling.  These  con¬ 
siderations  go  along  with  those  for  general  matrix  setup,  and  are  the  topic 
of  section  2.3.2.4. 

2.3.2A  Splitting  of  the  series  and  points  Recollecting  the  form  of  the  se¬ 
ries  in  section  2.3.1.2  (equation  31),  we  next  expand  the  pieces  of  equations 
(57)  through  (59)  and  arrange  the  terms  in  a  form  conducive  to  setting  up 
and  solving  the  linear  system  which  arises  from  collocation  of  a  boundary 
value  problem. 

Similarly,  the  collocation  points  in  figure  5  are  dealt  with  in  separate  parts. 
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to  facilitate  multiple  domain  handling  and  more  elaborate  unknown  repre¬ 
sentations^. 

For  collocation  point  ordering,  first  note  that  different  equations  are  valid  on 
each  of  the  boundaries,  and  further  that  not  all  of  the  ancillary  interior  points 
(triangles,  figure  5)  are  always  needed,  while  the  interior  points  (crosses,  fig¬ 
ure  5)  are  always  used.  This  suggests  separately  handling  the  boundaries, 
ancillary  interior  points,  and  interior  points,  giving  a  total  of  six  point  sets 
per  domain.  Throughout,  tiiese  will  be  denoted  by  t,b,l,r,a,i. 

For  the  unknown's  component  splitting,  note  the  following. 

•  The  Spline-Spline  product  terms  are  nonzero  in  most  of  the  point  re¬ 
gions. 

•  The  interior  series  has  zero  value  and  derivative,  or  zero  value  and  un¬ 
bounded  derivative  (depending  on  nuUifier)  on  the  boundaries,  and 
thus  never  needs  to  be  explicitly  evaluated  at  there. 

•  The  Series-Spline  terms  of  equation  (58)  and  (59)  are  nonzero  only  in 
the  interior  and  their  boimdary  (the  boimdary  on  which  the  spline's 
value  or  derivative  is  1) 

From  these  observations,  it  follows  that  for  multiple  domains,  only  the  Spline- 
Spline  product  terms  and  the  Series-Spline  terms  on  the  touching  boimd- 
aries  are  affected  by  the  overlap,  and  the  domains  unknowns'  are  otherwise 
independent.  Thus,  the  unknowns  are  split  into  top-value,  top-partial,  left- 
value,  . . .,  right-partial  and  interior  series  parts,  complemented  by  the  Spline- 
Spline  terms,  which  we  will  refer  to  as  comer  splines  (since  they  provide  val¬ 
ues  and  derivatives  at  the  comers). 

This  is  perhaps  best  illustrated  pictorially;  figure  6  shows  the  stmcture  of  a 
single  unknown  using  formulas,  while  figure  7  shows  the  stmcture  as  it  is 
implemented  on  the  computer. 

The  above  stmcturing  for  points  and  unknowns  leads  directly  to  the  matrix 
block  stmcture.  Since  this  stmcture  can  only  be  drawn  for  specific  cases,  we 
refer  here  to  figure  11  which  is  part  of  the  example  shown  later. 


2.3.2.5  Numerical  Example- Multiple  Unknown/Domain problem  This 
is  one  of  many  test  problems;  it  is  far  from  the  hardest.  It  was  chosen  here  be¬ 
cause  it  illustrates  all  important  features  of  the  sine  method  applied  to  multiple- 
domain,  multiple-unknown  P.D.E.  problems,  while  being  conceptually  sim¬ 
ple. 

®Mixed  boimdary  conditions  on  one  boimdary,  with  Dirichlet  or  Neumann  conditions 
on  another,  and  coupling  between  unknowns  and  domains. 
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Figure  6:  The  unknown's  representation  hierarchy  (mathematical  expres¬ 
sions  shown). 


The  Geometry 

This  problem  requires  the  solution  of  two  sets  of  elliptic  partial  differential 
equations  on  two  rectangles  which  have  one  common  botmdary.  See  figure 
8  for  the  equations  and  their  location.  In  the  equations,  the  following  defi¬ 
nitions  are  used: 

Gxx  —  (A  -H  2  G) 

=  (^  +  2G)^w  + A— w  (61) 
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Figure  7:  The  unknown's  representation  hierarchy  (program  parts  shown). 

and  etc.  are  determined  by  substitution  of  the  exact  answer,  chosen  below, 
into  the  appropriate  equation. 

The  exact  answers 

The  exact  answers  here  are  chosen  as 


ul  =  (x  —  (Xb  —  x) 

(63) 

(64) 

u2  =  {y-  Pa)  {Vb  -y)  +  (xb  - 

(65) 

w2  =  {y-  Pa)  (pb  -y){x-  Xa)  {Xb  - 

(66) 

and  runs  are  made  for  various  values  of  Ipul,  lpwl,.rpu2,  rpw2.  Xa  and  pa 
denote  the  left  and  bottom  boundary  of  the  domain,  respectively,  while  Xb 
and  pb  denote  the  right  and  top  boundaries. 

Numerical  Parameter  Values 

For  the  above  equations,  the  parameter  values  shown  in  tables  1  were  used. 


Parameters  of  Unknowns'  Sum  Representation 

Having  the  unknowns'  representation  as  detailed  in  section  2.3.2.4,  it  remains 
to  provide  parameter  values;  the  chosen  parameters  are  shown  in  table  2. 

Collocation  Points  Used 

In  figure  9,  the  used  collocation  points'  indices  are  shown  graphically.  Recall 
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O^xy  —  9l 
^yy  —92 


(^xy  —  p5 
^yy  ”  56 


^yy  —  53 
(^xy  —  54 


^xy  —  p7 

(Tyy  =  </8 


Figure  8:  Domains  with  equations. 
Table  1:  Equation  parameters 


domain  1: 

domain  2: 

Xa  =  1.0 

Xb  =  2.0 

Va  =  1.0 

Vb  =  2.0 

(9  =  1.1 

A  =  2.0 

Xa  =  2.0 

Xb  =  3.0 

Va  =  1-0 

Vb  =  2.0 
(9  =  1.1 

A  =  2.0 

that  the  actual  position  of  a  point  is  given  by 

Vk  =  (j)~^{kh2) 

and  interior  points  are  "bimched  up"  near  the  boundaries. 

The  block  matrix 

The  above  formulas  are  used  to  set  up  a  matrix  corresponding  to  the  origi¬ 
nal  linear  PDEs;  its  structure  is  shown  with  fixed  block  sizes  (to  keep  names 
legible)  in  figure  10,  and  in  proportion  in  figure  11 

Some  numerical  comparisons 

There  are  many  graphs  for  even  this  single  problem.  Some  representative 
results  are  shown  in  figures  13, 12  and  14.These  graphs  show  slices  in  the  j/ 
direction,  displayed  to  fit  on  single  pages.  Results  shown  use  29  series  terms 
in  both  X  and  y  directions,  and  were  computed  on  a  PC  with  32  Megabytes 


=  911 
=  S12 
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Table  2:  Parameters  for  the  sine  sum,  both  domains. 


A  I  I  « 


-Ml  Ni  -Ml  Ni 

Figure  9:  All  possible  collocation  points  are  shown,  but  only  boldface  ones 
are  used  in  this  example. 

of  memory. 

Because  the  accuracy  at  most  points  exceeds  2  digits,  most  graphs  show  no 
visible  errors;  the  partial  derivative  in  x  shown  in  figure  13  has  roughly  a 
10%  relative  error  on  the  y-slice  x  =  1.001  near  the  left  boxmdary  of  the  do¬ 
main.  Not  shown  is  the  relative  error  on  the  y-slice  at  a;  =  1 .003,  only  aroimd 
1%.  Put  another  way,  we  can  expect  less  than  1%  relative  error  when  within 
3/1000  of  the  singularity  or  crack,  using  sufficiently  many  terms  in  the  se¬ 
ries.  This  is  close  enough,  for  example,  to  curve  fit  a  simplified  local  model 
of  singular  behavior  near  the  crack  tip,  if  so  desired. 

In  conclusion,  we  thus  get  very  good  uniform  accuracy,  for  functions  and 
derivatives,  with  both  bounded  and  imboimded  derivatives,  making  the  method 
well  suited  for  crack  and  related  problems. 
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Figure  11:  Coupled  problem,  to  scale,  M  =  N  —  10. 
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Figure  12:  Function  value  comparison;  stars  are  computed  values,  solid  lines 
are  the  exact  answer. 


Sum/Exact  dCHiutin:  1  unknown:  u  Ml  «  14 


n - 

■■■1 

I^HII 

IIHI 

!■■■ 

I■■ll 

\mm\ 

■■■11 

immii 

1.001  1.111  1.222  1.333  1.443  1.5S4  1.665  1.775  1.886  1.997 

y  slices  at  indicated  x  values,  y  ranges  from  1.001  to  1.997 


Figure  13:  Comparison  of  function's  partial  derivative  in  x;  stars  are  com¬ 
puted  values,  solid  lines  are  the  exact  answer. 


Sum/Exact.  domain:  1  unknown:  u,  x  partial.  Ml  « 14 
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X  .  In  this  study,  the  loatbtransfer  characteristics  of  a  broken  fiber  are  investigated. 
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TT  Residual  stresses  due  to  curing  and  thermal  stresses  due  to  differences  between  the 
thermal  expansion  coefficients  of  the  matrix  and  fiber  may  have  a  major  effect  on  the 
microstresses  within  a  composite  material  system  and  must  be  added  to  the  stresses  induced  by 
the  external  mechanical  loads.  Such  microstresses  are  often  sufficient  to  produce  microcracking 
even  in  the  absence  of  external  mechanical  loads,  example  during  the  cooling  process. 


In  this  report  a  few  selected  results  are  presented  for  a  material  system  consisting  of  SIC-6 
cylindrical  fibers  which  are  periodically  embedded  into  a  plate  matrix  consisting  of  beta21 
material.  The  results  are  based  on  a  linear  elasticmicromechanics  model  which  provides  the  stress 
profiles  due  to  (i)  a  uniform  load  perpendicular  to  the  direction  of  the  fibers  and  (ii)  due  to  a 
thermal  expansion  mismatch. 


